A Computational Biomarker of Photosensitive Epilepsy from Interictal EEG

Abstract People with photosensitive epilepsy (PSE) are prone to seizures elicited by visual stimuli. The possibility of inducing epileptiform activity in a reliable way makes PSE a useful model to understand epilepsy, with potential applications for the development of new diagnostic methods and new treatments for epilepsy. A relationship has been demonstrated between PSE and both occipital and more widespread cortical hyperexcitability using various types of stimulation. Here we aimed to test whether hyperexcitability could be inferred from resting interictal electroencephalographic (EEG) data without stimulation. We considered a cohort of 46 individuals with idiopathic generalized epilepsy who underwent EEG during intermittent photic stimulation: 26 had a photoparoxysmal response (PPR), the PPR group, and 20 did not, the non-PPR group. For each individual, we computed functional networks from the resting EEG data before stimulation. We then placed a computer model of ictogenicity into the networks and simulated the propensity of the network to generate seizures in silico [the brain network ictogenicity (BNI)]. Furthermore, we computed the node ictogenicity (NI), a measure of how much each brain region contributes to the overall ictogenic propensity. We used the BNI and NI as proxies for testing widespread and occipital hyperexcitability, respectively. We found that the BNI was not higher in the PPR group relative to the non-PPR group. However, we observed that the (right) occipital NI was significantly higher in the PPR group relative to the non-PPR group. Other regions did not have significant differences in NI values between groups.


Introduction
People with epilepsy are prone to unprovoked epileptic seizures (Fisher et al., 2005). A subset of individuals also has seizures triggered by stimuli. The most common reflex epilepsy is photosensitive epilepsy (PSE), which occurs in ;2-5% of people with epilepsy (Padmanaban et al., 2019). In addition to seizures, visual stimuli can also induce EEG epileptiform activity in PSE. Thus, PSE can be diagnosed by recording a photoparoxysmal response (PPR) in an electroencephalogram (EEG) during intermittent photic stimulation (IPS). PSE is commonly seen in the idiopathic generalized epilepsies (IGEs), in particular, juvenile myoclonic epilepsy (JME), ranging from 30% to 90% (Wolf and Goosses, 1986;Appleton et al., 2000). PSE can serve as a useful model to understand epilepsy (Padmanaban et al., 2019). The PPR phenotype has been used within clinical trials to test the efficacy of antiseizure medication (Yuen and Sims, 2014) and in assessing other neurophysiological responses to visual stimuli (Parra et al., 2003;Perry et al., 2014). Thus, a better understanding of the pathophysiology of PSE may have an impact on the diagnosis and treatments of epilepsy.
Various studies have found evidence for both occipital and more widespread cortical hyperexcitability in people with PSE (Padmanaban et al., 2019). By using a pairedpulse flash-evoked technique, Strigaro et al. (2012) showed decreased inhibition in the visual system of people with PSE. Siniatchkin et al. (2007) used transcranial magnetic stimulation (TMS) to investigate the excitability of the visual and primary motor cortices in people with and without PPR. They observed that people with PPR have a significantly lower phosphene threshold (i.e., a lower TMS amplitude is required to elicit a perception of light in people with PPR). They also found that occipital TMS can more easily suppress visual perception in people with PPR than in those without it (Siniatchkin et al., 2007). Brigo et al. (2013) have also used TMS in PSE and provided further supporting evidence for these observations. Beyond the occipital lobe, there is evidence of the involvement of motor areas. Moeller et al. (2009Moeller et al. ( , 2013 have used EEG combined with functional magnetic resonance imaging (fMRI) to show increased blood oxygenation level-dependent (BOLD) signal within the visual cortex, and the premotor and parietal cortices after a PPR has been provoked by IPS. By using magnetoencephalography, Parra et al. (2003) observed enhanced network-level synchrony within the gamma band just before PPR formation. Vaudano et al. (2017) considered resting-state EEG-fMRI and found a decreased alpha-related inhibition of the visual cortex and sensory-motor networks in people with PSE versus other individuals with epilepsy (i.e., a lower anticorrelation between EEG alpha power and the BOLD signal in those networks). Diffusion tensor imaging studies have also revealed structural connectivity abnormalities within motor areas and the occipital lobe in individuals with PSE in IGE and JME, which may underpin hyperexcitability in those regions Vollmar et al., 2012).
Since a propensity for ictogenicity can be a consequence of hyperexcitability, here we aimed to find out whether widespread and/or occipital increases in the propensity for ictogenicity may be identified from resting interictal EEGs in people with PSE. As a proxy to evaluate network-wide and local propensity for ictogenicity, we used the brain network ictogenicity (BNI) and node ictogenicity (NI), respectively. BNI is a measure of how likely a functional brain network is to generate seizure-like activity in computer simulations (Petkov et al., 2014;Lopes et al., 2017). These simulations consist in placing a mathematical model of normal and abnormal epileptic EEG features into the functional network and computing the resulting brain dynamics. Functional networks more prone to seizures are expected to produce more seizure-like activity in the simulations (Petkov et al., 2014;Lopes et al., 2017Lopes et al., , 2021. Node ictogenicity is assessed by removing regions from the functional network and evaluating the resulting altered BNI Lopes et al., 2017). Brain regions whose removal produces a higher reduction of BNI are considered more likely to drive seizures.
We studied retrospective EEG recordings from two groups of individuals with IGE, one that had PPR during IPS (the PPR group), and another that did not have PPR (the non-PPR group). We tested the following two hypotheses: (1) the PPR group has a higher BNI than the non-PPR group; and (2) the PPR group has a higher occipital NI than the non-PPR group. We aimed to test whether interictal EEG may be used to predict PSE and whether the mechanisms of PPR are enduring features present in interictal brain states even without stimulation.

Data
We searched EEG reports at the University Hospital of Wales (Cardiff, UK) from 2007 to 2017 using the terms, "photoparoxysmal response," "PPR," "IGE," "JME," "JAE" (juvenile absence epilepsy), and childhood absence epilepsy "CAE," and limited our search to individuals who were 12-32 years of age at the time of EEG. Clinical EEG reports were reviewed to identify two cohorts of individuals with (1) IGE with IPS and PPR, and (2) IGE with IPS and no PPR. The IPS protocol at our center followed the international standard recommendations set out in the International League Against Epilepsy guidelines (Kasteleijn-Nolst Trenité et al., 2012). This comprised 10 s of IPS, 5 s with eyes open and 5 s with eyes closed, at the following incrementally increasing and then decreasing frequencies (1,2,4,6,8,10,12,14,16,18,20,60,50,40,30,25,and 20 Hz). The IPS was immediately stopped if a PPR was observed (i.e., a spike-wave discharge induced by the light stimulation). In some individuals, the order of stimulation was varied to look for reproducible PPR or to cease IPS if frequent PPR was already seen. We selected 29 individuals with PPR seen on EEGs and 20 individuals in the group without PPR. Henceforth, we will refer to the two groups as the PPR group and the non-PPR group, respectively.
The routine clinical EEG included ;20 min of restingstate interictal EEG preceding the IPS procedure. To test whether we could predict the individuals' responses based on the interictal EEG recorded before IPS, we selected three continuous segments of 20 s per individual within the 20 min window. The selection criteria were such as to avoid distinct (eye movement, muscle) artifacts and epileptiform activity, and such that all three segments were at least 1 min apart from each other. Three individuals in the PPR group were excluded from this study because we could not find interictal segments without marked artifacts within their EEG sessions. Thus, we studied a total of 46 individuals (26 PPR individuals and 20 non-PPR individuals). Tables 1 and 2 show the demographics and clinical information of these two groups. Figure 1 provides an example of a 20 s segment.
The EEG data were recorded at a sampling rate, f s , of 500 or 512 Hz. For consistency, we downsampled the data to 250 Hz. Furthermore, we applied a bandpass filter in the low alpha frequency band, 6-9 Hz (fourth-order Butterworth filter with forward and backward filtering to minimize phase distortions). We focused on the low alpha band because previous studies have shown that functional networks inferred from this frequency band are informative for both epilepsy diagnosis and epilepsy classification (Schmidt et al., 2016;Lopes et al., 2019), as well as of reduced inhibition in PSE (Vaudano et al., 2017). The data were also rereferenced to the average of all artifact-free segments.

Functional networks
To build functional networks from the 20 s interictal EEG segments, we followed a method based on previous studies (Schmidt et al., 2016;Lopes et al., 2019Lopes et al., , 2020Tait et al., 2021). We computed functional networks in the sensor space, where nodes corresponded to EEG channels. Connections between pairs of nodes i and j were inferred using the phase-locking value (PLV; Tass et al., 1998;Mormann et al., 2000), as follows: where N s is the number of time points (N t ¼ 5000), and Df ij ðt k Þ is the instantaneous phase difference between the EEG signals from channels i and j at time t k , which was computed using the Hilbert transform. We also calculated the average phase lag t ij between pairs of signals, as follows: We considered nodes i and j to be connected at PLV ij .0 and t ij .2p =f s with connection weight PLV ij . We neglected zero phase lag PLV (i.e., t ij ,2p =f s ) to avoid possibly artifactual relations because of volume conduction (Bastos and Schoffelen, 2016). Spurious connections were also rejected by comparing each possible PLV value to a set of PLV values obtained from surrogate time series (i.e., randomized time series comparable to the original time series). This comparison aims to remove connections whose PLV is because of random associations between signals and because of the finite nature of the signals. We generated 99 surrogates from the original EEG signals using the iterative amplitude-adjusted Fourier transform with 10 iterations Schmitz, 1996, 2000) and calculated 99 PLV values for every pair of channels. PLV values that are .95% of the corresponding PLV Figure 1. Summary of our computational method. We used interictal scalp EEG to infer functional networks. Then, to interrogate the networks, we considered a computational model of ictogenicity to simulate seizure-like activity on the functional networks. This allowed us to assess the propensity of a functional network to generate seizures in silico, the BNI. It further enabled us to measure the NI (i.e., the effect of removing a node on the BNI), which represents the local ictogenic propensity. values calculated from the surrogates were taken forward, representing the weights of the functional network. Thus, we constructed three functional networks per individual (i.e., one from each 20 s segment), each of them a matrix of statistically significant PLV values. Figure 1 shows an example of a functional network.

Computational model
We used these functional networks and a computational model of epilepsy to address two questions. (1) Do the networks from the PPR group have a higher ictogenic propensity than those from the non-PPR group? (2) Does the PPR group have a higher local ictogenic propensity in the occipital lobe compared with the non-PPR group?
The computational model was used to simulate the ability of a functional network to generate seizures, thereby enabling us to assess the propensity of the brain for ictogenicity and to compare it across individuals. It also enables us to analyze the relative importance of each brain region to the overall propensity to ictogenicity (i.e., the local ictogenic propensity; Lopes et al., 2017Lopes et al., , 2019Lopes et al., , 2020. Our approach consisted in using a phenomenological model of ictogenicity, the theta model (Ermentrout and Kopell, 1986), together with the calculated functional networks, to simulate brain-like activity and seizure-like transitions in these networks (Lopes et al., 2017(Lopes et al., , 2019(Lopes et al., , 2020. In this model, the activity at network node i is represented by a phase oscillator, u i , that obeys the following ordinary differential equation, as follows: where I i t ð Þ is an input current to the node i at time t. This current is given by the following: where I 0 1j i ð Þ t ð Þ is Gaussian noise, K is a global scaling factor of network interactions, and N is the number of nodes in the network (N ¼ 20). The noise represents signals coming from brain regions outside of the considered functional networks and the term ½1 À cosðu j À u s ð Þ Þ is the output of node j, representing a displacement from its resting-state phase u s ð Þ . PLV ji is j, the ith entry of the adjacency matrix that represents the functional network calculated above. Within the model, oscillators may transit between the following two states: a resting state characterized by fluctuations close to the stable phase u ðsÞ if I i ,0, and a seizure-like state represented by a rotating phase if I i .0. The transition between the two states corresponds to a saddle node on an invariant circle bifurcation at I i ¼ 0. This oscillator model was shown to approximate a neural mass model for the purpose of studying the ictogenicity of a network (Lopes et al., 2017). Figure 1 shows an example of simulated signals using this model. We chose parameters according to previous studies (Lopes et al., 2017(Lopes et al., , 2019(Lopes et al., , 2020: I 0 ¼ À1:2 and noise SD s ¼ 0:6. The global scaling factor K was treated as a free parameter (see the Brain network ictogenicity and node ictogenicity section).

Brain network ictogenicity and node ictogenicity
We measured the ability of the brain network to generate seizures by using the concept of BNI. The BNI is the average fraction of computational time that the network nodes spend in the seizure-like state (Petkov et al., 2014;Lopes et al., 2017), as follows: where t i ð Þ sz is the time that node i spends in the rotating state during a total simulation time T. We used T ¼ 4 Â 10 6 , as in previous studies (Lopes et al., 2019(Lopes et al., , 2020; but see Lopes et al., 2017, for more details about the calculation of t i ð Þ sz ). To avoid an arbitrary choice of the free parameter K, we considered a robust redefinition of the BNI (Lopes et al., 2017 given by the following: where K 1 and K 2 were chosen so that to capture the full variation of the BNI from 0 to 1 for all networks under consideration. We used the same interval K 1 ; K 2 ½ ¼ ½1; 40 for all functional networks of all individuals. For each individual, we obtained three d BNI values, one per each functional network, and took their mean value for the statistical analysis below.
To quantify the relative importance of each node to the ability of the network to generate seizures (i.e., the local ictogenicity, we considered the NI; Goodfellow et al., 2016;Lopes et al., 2017Lopes et al., , 2019. The NI i ð Þ measures the effect of removing node i on the ability of the network to generate seizures in silico, and it is given by the following: where BNI pre is the BNI before node removal, and BNI i ð Þ post is the BNI after the removal of node i. Following previous studies, we selected the parameter K such that BNI pre ¼ 0:5 Lopes et al., 2017Lopes et al., , 2019Lopes et al., , 2020. BNI i ð Þ post is typically smaller than BNI pre , meaning that NI i ð Þ is usually positive. Thus, the NI i ð Þ ranges typically between 0 and 1, where 0 corresponds to a node removal that has no effect on seizure generation (BNI i ð Þ post ¼ BNI pre ), and 1 corresponds to a node removal that stops seizures in the network (BNI i ð Þ post ¼ 0). The higher the NI, the more important the node is for seizure dynamics. As above, we also took the mean NI values of each node across the three functional networks per individual for the statistical analysis presented below. Figure 1 summarizes the key steps of our methods.

Statistical methods
We used a one-sided Mann-Whitney U test to assess whether the mean BNI values were higher in the PPR Research Article: New Research group than in the non-PPR group. We also used the same test to compare the mean NI in the occipital areas of the two groups. Furthermore, we performed an exploratory analysis where we used the two-sided Mann-Whitney U test to evaluate whether the mean NI of each brain area (other than the occipital areas) was different between the two groups. To correct for multiple comparisons within each family of tests, we applied the Bonferroni-Holm procedure to each of these three analyses separately. For all of these Mann-Whitney U tests, we report the U-statistic and z score. We further quantified statistical differences by using estimation statistics (Ho et al., 2019). We used the median difference as the effect size and 95% confidence intervals (CIs) were built using 5000 bootstrap samples. To measure these estimation statistics, we used the tools provided by (Ho et al., 2019; https://www.estimationstats.com).
In the case of statistically significant results, we further quantified the difference between the two groups using the receiver operating characteristic (ROC) curve, area under the curve (AUC), and sensitivity and specificity.

Data availability
MATLAB scripts implementing the methods described in the article are freely available online at https://github. com/ml0pe5/Photostimulation_BNI_NI.

Results
Do people with epilepsy who had a PPR to IPS have higher widespread excitability than people who did not have a PPR to IPS? Assuming that they do have a higher widespread excitability, does it translate into a higher propensity to seizures as assessed from their resting-state functional networks? To address this latter question, we measured the BNI of each individual from interictal EEG recordings before the IPS study. Figure 2 shows the BNI of the two groups: the PPR group (i.e., positive PPR on IPS) and the non-PPR group (i.e., no PPR on IPS). We found that the BNI is not statistically significantly higher in the PPR group (p ¼ 0:89; one-sided Mann-Whitney U test, U ¼ 556, z ¼ À1:23). We further found that the unpaired median difference of BNI was À1:48 with a 95% confidence interval equal to -4.64 to 1.40. These results suggest that PPRs are not a consequence of an individual having a higher propensity to generate seizures (i.e., a lower threshold for seizure emergence).
Having observed that the PPR group is not more prone to seizures than the non-PPR group, we hypothesized that their paroxysmal response to light stimulation results from a higher ictogenic propensity in the occipital lobe compared with the non-PPR group. Thus, we tested whether the occipital areas of PPR individuals had a stronger contribution to seizure generation within the whole brain network than they did in non-PPR individuals. Figure 3 compares the NI in the occipital areas (channels O1 and O2) between the PPR and non-PPR groups. We observed that the median of the NI of the O1 region was higher in the PPR group relative to the non-PPR group, but it was not significant (p ¼ 0:18; one-sided Mann-Whitney U test, corrected for multiple comparisons, U ¼ 652, z ¼ 0:90). We found a significantly higher NI in the O2 region in the PPR group relative to the non-PPR group (p ¼ 0:013; one-sided Mann-Whitney U test, corrected for multiple comparisons, U ¼ 724, z ¼ 2:49). Estimation statistics further supported these observations: for O1, the unpaired median difference of NI was 0.03 with a 95% confidence interval of À0.02 to 0.06; and for O2 the effect size was 0.07 with a 95% confidence interval of 0.01-0.12. As an exploratory analysis, we also tested whether any other region had a different ictogenicity in the PPR group compared with the non-PPR group, and we found that no region had a statistically different NI between the two groups (Fig. 4, Table 3).
Given the observation that the O2 was the most different area in terms of NI between the two groups, we further measured its ability to classify the two groups. We calculated the ROC curve and found an AUC of 0.72, a sensitivity of 0.85, and a specificity of 0.55.
To better characterize our findings, we performed a supplementary analysis. Considering that previous studies have suggested that photosensitive epilepsy may be differentiated from other types of epilepsy by analyzing resting-state EEG power spectrum (Vaudano et al., 2017) or by analyzing resting-state EEG functional networks (Varotto et al., 2012), we tested whether such types of analyses could distinguish our two groups of individuals. Our purpose was to determine whether our findings could be explained by these straightforward spectral or connectivity analyses of EEG signals. Based on the study by Vaudano et al. (2017), we computed the relative power in the alpha band in the occipital electrodes (O1 and O2) and tested whether the relative power was higher in the PPR group than in the non-PPR group. We observed no statistical differences (Table 4). Furthermore, we tested whether the occipital electrodes had higher connectivity strength in the PPR group relative to the non-PPR group (Rubinov and Sporns, 2010). Again, we found no statistical differences (Table 4).

Discussion
In this study, we consider the hypothesis that PSE may be underpinned by both occipital and more widespread cortical hyperexcitability (Padmanaban et al., 2019). We assess whether such hyperexcitability is an enduring feature inferable from resting brain activity of people with PSE. We used interictal scalp EEG recordings from two groups of individuals with IGE, one that had PPR during IPS, and another that did not have PPR. We inferred functional networks from the EEG data and studied the propensity of the participants' brain networks to generate seizure activity using a computer model. By simulating the emergence of seizure-like dynamics on the networks, we were able to compare the overall ictogenic propensity and local ictogenic propensity across individuals.
We first tested whether the PPR group had a higher BNI than the non-PPR group. A higher BNI would suggest a higher propensity of the PPR group to have seizures, which could be indicative of widespread hyperexcitability in this group. We found that the BNI was not statistically different between the two groups. While this result suggests that there is not an enduring hyperexcitability identifiable from interictal EEG data, we acknowledge that there is an alternative interpretation. Namely, the BNI may better reflect the expectable seizure frequency of an individual  rather than their hypothesized widespread hyperexcitability. From this perspective, the BNI was not expected to be higher in the PPR group relative to the non-PPR group, as observed, because people with PSE do not have higher seizure rates than people without PSE (Covanis, 2005;Verrotti et al., 2012).   Fig. 3). p-Values, effect sizes, and confidence intervals are presented in Table 3. All these differences in NI are not statistically significant.
We then tested whether the PPR group had a higher occipital NI than the non-PPR group. We found that indeed the NI is higher in occipital regions of the PPR group relative to the non-PPR group, in particular in the right occipital region (O2). This result suggests that there is an enduring hyperexcitability of the occipital regions in PSE, and that this hyperexcitability can be inferred even from apparently normal scalp EEG recordings, without requiring the use of stimulation. More specifically, a higher NI of a region indicates that this region is more relevant to seizure emergence. Note, however, that we do find other regions with higher NI than the occipital lobe (e.g., the frontal areas, FP1 and FP2; Fig. 4). Together, these results suggest that people with PSE have a higher ictogenic propensity in their occipital lobe compared with other individuals with epilepsy, but that this region is not their main driver of seizures (in line with people without PSE). We further showed that, using NI to classify people with PSE, we obtained a sensitivity of 0.85 and a specificity of 0.55. Our results are in line with previous evidence that has indicated occipital hyperexcitability in people with PSE (Wilkins, 1995;Porciatti et al., 2000;Ferlazzo et al., 2005;Padmanaban et al., 2019).
There have been previous studies assessing PSE based on resting-state EEG functional connectivity. Varotto et al. (2012) computed partial directed coherence from 10 individuals with PSE and 10 healthy control subjects. They found that people with PSE had enhanced connectivity, predominantly involving the anterior cortical regions. A key advantage of our study compared with this one is that we compared two groups of people with epilepsy, one with PSE and one without, allowing us to assess features that should be specific to PSE, rather than of epilepsy more generally. Vaudano et al. (2017) analyzed resting-state EEG-fMRI from 16 individuals with genetic generalized epilepsy with PSE, 13 individuals with genetic generalized epilepsy without PSE, and 15 individuals with focal epilepsy. They found that the PSE group had significantly higher mean alpha power than the other two groups. Based on both EEG and fMRI, they showed that the "cortical-subcortical network generating the alpha oscillation at rest is different in people with epilepsy and visual sensitivity." While not having fMRI in our dataset, our analysis went The relative power was computed following the methods of Vaudano et al. (2017). We considered the following two frequency bands: low alpha power (6-9 Hz) as in the main analysis; and alpha power (7.5-12.5 Hz) as in the study by Vaudano et al. (2017). The occipital connectivity corresponds to the connection strength of the electrodes (i.e., sum of in-strength and out-strength; Rubinov and Sporns, 2010). The p-values, U statistics, and z-scores correspond to one-sided Mann-Whitney U tests assessing whether the relative power (or connectivity strength) is higher in the PPR group relative to the non-PPR group at a given occipital electrode. All p-values are not significant. The effect size (median difference) and confidence intervals further show that there are no statistical differences between the groups when using these measures. The p-values, U statistics, and z-scores correspond to two-sided Mann-Whitney U tests assessing whether NI is different between the two groups at a given node. The p-values were corrected using the Bonferroni-Holm procedure, and none is significant. The effect size (median difference of NI) and confidence intervals further show that there are no statistical differences between the groups in any of these nodes.
beyond those of these two studies. Rather than focusing on differences in network structure between groups, we placed a model of ictogenicity on the networks to pose hypotheses about how their structure may influence ictogenicity. Additionally, we tested whether we could distinguish our two groups of individuals based on analyses of alpha power and connectivity strength, similar (but not equivalent) to those performed by Vaudano et al. (2017) and Varotto et al. (2012). We found that those measures could not differentiate our groups of individuals, which further supports the use of our computational method. Nevertheless, we note that our analyses of EEG spectrum and connectivity are not replications of the analyses reported by Vaudano et al. (2017) and Varotto et al. (2012). In addition to some methodological differences, both the data and cohorts of individuals were not directly comparable to ours. Vaudano et al. (2017) used 32-channel EEG and fMRI, while we used 20-channel EEG. Varotto et al. (2012) compared individuals with PSE with healthy control subjects, whereas we assessed a cohort of IGE patients with and without PSE. Our purpose was not to replicate their studies, but rather to assess whether we could differentiate our two groups of individuals using straightforward EEG spectrum or connectivity analyses, which to some extent underlie our computational framework. Future work should aim to replicate together our study and their studies to better assess whether the results are concordant or complementary.
In conclusion, we applied the BNI and NI frameworks to interictal scalp EEGs from people with epilepsy, and we found that people with PSE have higher NI in occipital regions than people without PSE. Our findings suggest that the mechanisms of PSE may be revealed by enduring features present in interictal brain states. Our results also suggest that the NI applied to resting-state EEG may aid the diagnosis of PSE without the need of stimulation.